package de.hhn.mim.psb;

/**
 * 
 */

import java.util.Arrays;

import Jama.Matrix;
import Jama.EigenvalueDecomposition;
import Jama.SingularValueDecomposition;

/**
 * Rumpf-Klasse zur Implementierung einer Hauptkomponentenanalyse
 * Grundlage der Hauptkomponentenanalyse ist die Aufstellung der
 * Kovarianzmatrix, deren Eigenwerte und Eigenvektoren bestimmt
 * werden.
 * Die Bestimmung der Eigenwerte und Vektoren erfolgt mit Hilfe von
 * Methoden aus der Jama Bibliothek (s.u.) in der Methode 
 * calcEigenvalueDecomposition. Die Berechnung der Kovarianzmatrix
 * muss zuvor in der Methode calcCovMatrix erfolgen.
 * Diese Methode muss implementiert werden. Sie gibt zurzeit keine
 * Matrix zurueck, als Ergebnis stuerzt das Programm ab.
 * 
 * berechnet Eigenwerte und Eigenvektoren auf Basis einer 
 * Hauptkomponentenanalyse fuer eine Menge gegebener Punkte vom Typ Vect3D
 * Klasse kapselt Funktionalitaet aus Jama Bibliothek
 * http://math.nist.gov/javanumerics/jama/
 * 
 * Die Berechnung der Eigenwerte und Vektoren erfolgt beim Aufruf des 
 * oeffentlichen Konstruktors, auf Basis der �bergebenen Punkte
 * Das Ergebnis kann mit getEigenValues und getEigenVectors abgerufen werden
 * getSortEigenVals und getSortEigenVects liefern die Eigenwerte und Vektoren
 * sortiert nach dem groessten Eigenwert zurueck.
 * calcCovMatrix liefert die Kovarianzmatrix, fuer einen gegebenen Array von 
 * Punkten (Vect3D) 
 * calcCenterOfMass liefert den Schwerpunkt fuer die Punkte eines gegebenen
 * Arrays.
 * 
 * @author Rolf Bendl
 *
 * 2012-05-12 Klasse erstellt
 */
public class PCA3D
{
	Vect3D []   points;
	double [][] covMatrix;
	double [][] eigenVectors;
	double []   eigenValues;
	double []   sortEigenVals;
	double [][] sortEigenVects;
	/**
	 * 
	 */
	private PCA3D()
	{
		points           = new Vect3D[0];
		covMatrix        = new double[3][3];
		covMatrix[0]     = new double[0];
		covMatrix[1]     = new double[0];
		covMatrix[2]     = new double[0];
		eigenVectors     = new double[3][3];
		eigenValues      = new double[3];
		sortEigenVects   = new double[3][3];
		sortEigenVals    = new double[3];
	}

	/**
	 * constructor
	 * @param points zur PCA Berechnung
	 */
	public PCA3D(Vect3D [] points)
	{
		this();    
		this.points = new Vect3D[points.length];
		for(int i = 0; i < points.length; i++)
			this.points[i] = new Vect3D(points[i]);

		covMatrix = calcCovMatrix(points);
		calcEigenvalueDecompostion();
		sortEigenValsAndVectors();
	}
	/**
	 * @return die Eigenwerte
	 */ 
	double [] getEigenValues()
	{
		double [] result = new double [eigenValues.length];
		for(int i = 0; i < eigenValues.length; i++) result[i] = eigenValues[i];
		return result;
	}
	/**
	 * @return die Eigenvektoren
	 */
	double [][] getEigenVectors()
	{
		double[][] result = new double [eigenVectors.length][eigenVectors[0].length];
		for(int i = 0; i < eigenVectors.length; i++)
			for(int j = 0; j < eigenVectors[0].length; j++)
				result[i][j] = eigenVectors[i][j];
		return result;
	}

	/**
	 * liefert Array mit den Indices der Eigenwerte in absteigender 
	 * Ordnung (groesster zuerst)
	 * @return Array mit Indices
	 */
	private int [] calcEigenValuesOrder()
	{
		double [] vals = getEigenValues();
		int [] result = new int [vals.length];
		for(int i = 0; i < vals.length; i++)
		{
			double largeV = Double.MIN_VALUE;
			int    largeJ = 0;
			for(int j = 0; j < vals.length; j++)
			{
				if(vals[j] > largeV)
				{
					largeV = vals[j];
					largeJ = j;
				}
			}
			result[i]    = largeJ;
			vals[largeJ] = Double.MIN_VALUE;        
		}
		//    for(int i = 0; i < vals.length; i++)
			//      System.out.println(result[i]);

		return result;
	}

	/**
	 * sortiert die Eigenwerte und Vektoren und speichert
	 * sie in den Variablen sortEigenVals und sortEigenVects
	 */
	private void sortEigenValsAndVectors()
	{
		int []      order = calcEigenValuesOrder();
		double []   vals  = getEigenValues();
		double [][] vects = getEigenVectors();

		for(int i = 0; i < vals.length; i++)
		{
			int idx = order[i];
			sortEigenVals[i] = eigenValues[idx];
			for(int j = 0; j < eigenVectors[idx].length; j++)
			{
				sortEigenVects[j][i] = eigenVectors[j][idx];
			}
		}
	}
	/**
	 * @return die sortierten Eigenwerte (groesster zuerst)
	 */
	double [] getSortEigenVals()
	{
		double [] result = new double [sortEigenVals.length];
		for(int i = 0; i < sortEigenVals.length; i++) result[i] = sortEigenVals[i];
		return result;
	}
	/**
	 * @return die sortierten Eigenvektoren (Vektor zum groessten EW zuerst
	 */
	double [][] getSortEigenVects()
	{
		double[][] result = new double [sortEigenVects.length][sortEigenVects[0].length];
		for(int i = 0; i < sortEigenVects.length; i++)
			for(int j = 0; j < sortEigenVects[0].length; j++)
				result[i][j] = sortEigenVects[i][j];
		return result;
	}
	/**
	 * berechnet Eigenwerte und Eigenvektoren
	 */

	private void calcEigenvalueDecompostion()
	{
		Matrix                  jamaMat     = new Matrix(covMatrix);
		EigenvalueDecomposition evd         = jamaMat.eig();
		double []               eigenVals   = evd.getRealEigenvalues();
		Matrix                  eigenVects  = evd.getV();
		Matrix                  eigenValMat = evd.getD();
		int                     print       = 0;
		eigenVectors = eigenVects.getArrayCopy();
		eigenValues  = new double[3];
		for(int i = 0; i < eigenVals.length; i++) eigenValues[i] = eigenVals[i];

		if(print > 1)
		{
			for(int i = 0; i < covMatrix.length; i++)
			{
				System.out.printf("covMatrix[%d][] ", i);

				for(int j = 0; j < covMatrix[i].length; j++)
				{
					System.out.printf("%8.4f ", covMatrix[i][j]);
				}
				System.out.println("\n");
			}
		}
		if(print > 1)
		{
			System.out.println("jamaMat");
			jamaMat.print(8, 3);
		}
		if(print > 0)
		{
			System.out.println("eigenVects");
			eigenVects.print(8,3);
			System.out.println("eigenValMat");
			//      for(int i = 0; i < eigenVals.length; i++) 
			//        System.out.printf("%8.3f\n",eigenValues[i]);
			eigenValMat.print(8, 3);
		}
		// ergibt die gleichen Eigenwerte, aber in einer anderen Reihenfolge
		/*    
    SingularValueDecomposition svd      = jamaMat.svd();
    double []                  svals    = svd.getSingularValues();
    Matrix                     svds_mat = svd.getS();
    Matrix                     svd_umat = svd.getU();
    Matrix                     svd_vmat = svd.getV();

    System.out.println("SVD S-Matrix");
    svds_mat.print(8,3);
    for(int i = 0; i < svals.length; i++) 
       System.out.println(svals[i]);
    System.out.println("SVD U-Matrix");
    svd_umat.print(8,3);
    System.out.println("SVD V-Matrix");
    svd_vmat.print(8,3);

    // liefert bereits sortierte Werte
    eigenValMat  = svds_mat;
    for(int i = 0; i < svd_umat.getRowDimension(); i++)
      for(int j = 0; j < svd_umat.getColumnDimension(); j++)
        eigenVectors[i][j] = svd_umat.get(i, j);

    for(int i = 0; i < svals.length; i++) 
      eigenValues[i]  = svals[i];
		 */    
	}

	/**
	 * Berechnet Kovarianzmatrix
	 * 
	 * Lt. Script MH-03 Teil 1 Segmentierung mit Form- und Erscheinungsmodellen:
	 * gegeben ist die Datenmatrix durch Anordnung der Punkte in Spalten
	 * Matrix M ergibt sich durch Subtraktion des Mittelwertes von jeder Zeile
	 * 1/3 * M * M-transpo ist die Kovarianzmatrix
	 * @param points
	 */
	double [][] calcCovMatrix(Vect3D [] points)
	{
		int         length     = points.length;
		double [][] values     = new double[3][length];
		int         print      = 0;

		System.out.println("Achtung, die Methode PCA3D::calcCovMatrix ist nur " + 
				"eine Dummy Implementierung, \nSie muessen die " +
				"Berechnung der Kovarianzmatrix implementieren");

		// Achtung nur Dummy-Implementierung, damit die Methode eine
		// gueltige Matrix zurueck gibt, und das Programm nicht abstuerzt
		//
		// Hier muessen Sie die Berechnung der Kovarianzmatrix implementieren
		//
		// Die Kovarianzmatrix darf nicht vom Typ Matrix3D sein, dieser
		// Typ realisiert ausschliesslich eine 4x4 Matrix
		//
		// Wenn Sie allgemeine Matrix-Operationen benoetigen, verwenden Sie
		// den Typ Matrix aus der Jama-Bibliothek
		// 
		Vect3D m = calcCenterOfMass(points);
		
		
		System.out.println("mean: ("+m.x+","+m.y+","+m.z+")");
		//mittelwert von jeder Spalte abziehen	
		for(int i = 0; i < length; i++){
			values[0][i] = points[i].x-m.x;
			values[1][i] = points[i].y-m.y;
			values[2][i] = points[i].z-m.z;
		}
		Matrix M = new Matrix(values);
		M = M.times(M.transpose());
		
		M = M.times(1.0/points.length);	
		
		//umwandeln von Matrix in 2D-Array und Rueckgabe
		return M.getArray();
	}

	/**
	 * Berechnet den Schwerpunkt einer gegebenen Punktmenge
	 * @param points, die Punkte
	 * @return den Schwerpunkt
	 */
	public Vect3D calcCenterOfMass(Vect3D [] points)
	{
		Vect3D center = Vect3D.mean(points);
		return center;
	}
}
